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We consider spherically trapped Bose gases in three dimensions with contact interactions, and 
investigate whether the Bose-Einstein condensate at zero temperature is stable against macroscopic 
fragmentation into a small number of mutually incoherent pieces. Our results are expressed in 
terms of a dimensionless interaction measure proportional to the Thomas-Fermi parameter. It is 
shown that while three-dimensional condensates are inherently much more stable against macro- 
scopic fragmentation than their quasi-one- and quasi-two-dimensional counterparts, they fragment 
at a sufficiently large value of the dimensionless interaction measure, which we determine both fully 
numerically and semi-analytically from a continuum limit of large particle numbers. 
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I. INTRODUCTION 

Bose-Einstein condensation (BEC) 0, Q of noninter- 
acting bosons can in principle occur in arbitrarily large 
spatial dimension D, depending both on the properties of 
the single-particle spectrum and the confining potential 
Q. On the other hand, the Hohenberg-Mcrmin- Wagner 
theorem rules out, independent of the strength of interac- 
tions, BEC in D < 2 The latter theorem, however, 
applies to homogeneous condensates in the thermody- 
namic limit, where long-range phase fluctuations trigger 
the decay of the BEC into infinitely many fragments. 
This changes for trapped condensates, where the finite 
extension of the gas cuts off the phase fluctuations in the 
corresponding directions. While an explicitly interaction- 
independent formulation of the theorem is still possi- 
ble, the geometric shape of the condensate enters the 
Bogoliubov inequality on which the Hohcnbcrg-Mcrmin- 
Wagncr theorem rests [(|. 

It is well known that in three spatial dimensions, frag- 
mentation does not occur in the thermodynamic limit 
in a homogeneous system for positive coupling 0, H|, 
while for negative coupling constant the system is un- 
stable. For a trapped, that is spatially localized and 
inhomogencous system, in three dimensions, the rele- 
vant dimensionless parameter to measure the importance 
of interactions over the single-particle kinetic and trap- 
ping contributions to the energy is the Thomas-Fermi 
parameter Na s /lo (where a s and Iq are s-wave scatter- 
ing and harmonic trapping length, respectively). It was 
previously observed by us that harmonic trapping and 
positive interaction coupling can lead to fragmentation 
into two mutually incoherent macroscopic pieces, form- 
ing a so-called fragmented condensate, well before the 
thermodynamic limit is taken for quasi-one-dimcnsional 
(quasi-lD) and quasi-two-dimensional (quasi-2D) gases 
@, Q3 ■ To more completely elucidate the dimension de- 
pendence of the many-body physics of fragmentation, we 
present here a detailed analysis of fragmentation for the 



completely symmetric example a spherically trapped 3D 
condensate. Viewed from a different angle, we investi- 
gate to which extent the conventional textbook wisdom 
[111 [l2| , that when the Thomas- Fermi parameter of a 3D 
spherically trapped condensate is going to infinity yields 
a interaction-dominated single Bose-Einstein condensate 
(with parabolic shape in this Thomas-Fermi limit and in 
a harmonic trap) needs revision. 

In three dimensions, due to the spherical symmetry 
of the system, even when the field operator expansion 
is restricted to the low-energy sector, there are poten- 
tially four single-particle states which are macroscopi- 
cally occupied. By numerical analysis and general sym- 
metry arguments, we find that fragmentation is domi- 
nated by two orbitals at a dimensionless coupling mea- 
sure which is proportional to the Thomas-Fermi param- 
eter. The critical coupling measure is one (two) orders 
of magnitude larger than the corresponding measure in 
the quasi-2D (quasi-lD) cases. In addition, the maximal 
degree of fragmentation @ turns out to be significantly 
smaller than in the latter cases. Our result therefore im- 
plies the rapidly growing persistence of an interacting, 
trapped scalar Bose-Einstein condensate against macro- 
scopic fragmentation upon increasing the spatial dimen- 



II. SPHERICALLY TRAPPED GASES IN 
THREE DIMENSIONS 

A. The four-mode approximation for the 
Hamiltonian 

To facilitate comparison with the previously treated 
quasi-lD and quasi-2D trapping cases, we will make a 
one-parameter variational ansatz for the single-particle 
orbitals like in 0] , This involves ground and first excited 
states of the harmonic oscillator, with the variational pa- 
rameter chosen to be harmonic oscillator length. Com- 
pared to fully self-consistent multiconfigurational Hartrec 
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calculations as performed, e.g., in [13M16I ] . while be- 
ing less quantitatively accurate, the variational approach 
leads to a qualitatively correct picture of the fragmenta- 
tion phenomenon. A particular merit of this approach 
is that the parameter dependence of the fragmentation 
transition is transparent: We find that fragmentation is 
decided by a single parameter, G3 in Eq. below, which 
measures the relative importance of interactions over the 
single-particle (trapping) energies. In addition, the vari- 
ational approach is capable to deal with the limit of very 
large particle numbers ./V; in its continuum limit, which 
we will derive below, there is indeed no upper bound to 
the value of N. This is particularly beneficial in three 
spatial dimensions, where the particle numbers at the 
same densities are obviously larger than in one- and two- 
dimensional systems; we were able to numerically calcu- 
late within relatively short timcscalcs systems with up to 
N ~ 10 6 particles. 

To formulate the proper variational orbitals basis, we 
first write down the well-known eigenstates and energies 
of the isotropic harmonic oscillator in three spatial di- 
mensions, 
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with normalization N n 1 = ^-i"^' 1 ^ , and the functions 

' Vn!r(!+5) 

Y™, Ln +1 ^ 2 ^ are spherical harmonics and generalized La- 
guerre polynomials respectively. 

The radial distance r is scaled by a length R, which 
for noninteracting condensates is given by the harmonic 
oscillator length, R = l = uj^ 1 / 2 (h = M = 1, where M 
is the boson mass). In the following, we assume R to be a 
variational parameter, which will determine the family of 
solutions of the many-body equations, i.e., whether single 
or fragmented condensates are obtained is determined by 
variation of R. 

According to J]}, the four energetically lowest states 
are given by the quantum numbers n = 0, / = (ground) 
and n = 0, I = 1, m = —1,0,1 (first excited) [adopting 
the Condon-Shortlcy phase convention for spherical har- 
monics! , 
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The contact-interaction many-body Hamiltonian, ob- 
tained by truncating the field operator expansion after 



including the four modes ([2]), reads 
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Note that the pair-exchange scattering (terms cx Ei) oc- 
curs also between energetically degenerate orbitals (de- 
generate on the single-particle level), being represented 
by the term oc E3, and involving the excited states 
m = ±1 and m = (I = 1). This is distinct from the 
quasi-lD and quasi-2D cases treated in where pair- 
exchange scattering only occurs between pairs of ground 
and excited single-particle states. 

The interaction matrix elements Vijki = 
gj J J r 2 sin 6 drd6d(j)ipi (r)ip*(r)il>k(i')'4>i( r ) arc related 
to the coefficients in ([3]) as follows. The no n- vanishing 
pair-exchange matrix elements are 



Ei — V1100, 
Es = — 11 



E-2 = V+-ov + V- 
V-+u- 



-00, 



(4) 



The remaining coefficients are of the density-density 
type, 

Co = Voooo, Ci — Vim, C+ = 

C_ = V , D\ = V0101 + V1010 + V1001 + KniOi 

D 2 = D 1 (l->-±) ) D 3 = Ui(0 -> ±), 

D 4 = D 1 (0->-+ > l ->•-). (5) 

The result for the scattering coefficients (reinstating now 
the variational harmonic oscillator length) , may be writ- 
ten in compact vector notation in the following way, 



{Co, Ci, C±, Di, D 2 , D3, D4, Ei,E 2 , E3} 
9 
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The single-particle energies are given by ti 
J r 2 sin edrd9d(j) [\V^i\ 2 /2 + w 2 r 2 |^| 2 / 2 ] , and read 
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Defining the scaled variational parameter A = R/Iq, we 
have as the typical units of single-particle energies and 
coupling constants 
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where we introduced the dimensionless interaction cou- the single parameter deciding the question of coherence 
pling versus fragmentation. 



Ng 

G3 = — — vs/2) . (9) B. Wavefunction ansatz and the eigenvalue 

problem 



(27T)3/2Z ' 



Like its quasi-lD and quasi-2D counterparts G\ = 

(dp7% and G2 = (2i!)z%i- ' wnere k and l± arc nar " 
monic oscillator lengths of a cylindrical trap, the quantity 
G3 measures the relative importance of total interaction 
and kinetic-potential energy terms in the Hamiltonian. 
In the present spherically trapped 3D case, G3 is simply 
directly proportional to the well-known Thomas-Fermi 
parameter [12j, which, as will be demonstrated below, is 



We employ a general four-mode ansatz in the Fock 
subspace of fixed total particle number 

|*>= Y,^M,i-\ N -h-l+-l-,hJ+, 1-). (10) 

The total energy E = (^f\H \^f) in terms of the level 
occupation amplitudes then reads 



E = c h,i± E l^i.'il 2 + l E i d i^h,i+,i-^h+2,i + ,i- + -E 2 d 2 <>Pi ul+tl _ipi 1 ,i + +i,i_+i + ^E 3 d 3 tpf ul+!l _'>pi 1 -2,i + +i,i_ + i 
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where the diagonal and pair-exchange coefficients take the explicit form 
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Finally, the minimization of the energy functional with respect to ^* l+ l gives the eigenequations 

E\ E 2 
Eipi u i + ,i_ = ci u i ± ipi u i + j_ + — d 1 (h,l + ,l-)il>i 1+ 2,i + ,i_ + — d2(li,l+,l-)ipi 1 ,i + +i,i_+i 

QjQ E\ 

+ -^-^3^1^+) l-mi-2,l + +l,U+l + —dl{h - 2, l+, l-)lpl x -2,l+,l- + —<k(li,l + -1,1-- lj^i.i-,--!,!--! 
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III. SOLVING THE EIGENVALUE PROBLEM cause the matrix dimensions scale roughly with N 3 x N 3 

when naively implemented. However, the assumed or- 



A. Decomposition into smaller problems: The 
fc-subspaces 

We are facing a high-dimensional eigenvalue problem 
that is difficult to solve for typical particle numbers be- 



bitals allow for algebraic simplifications, to be explained 
in what follows, in order to significantly reduce the prob- 
lem size. 

Observe that only a particular set of couplings be- 



4 



tween the l\ and l± terms appear in cq. (fT3|) . To 
be more specific, only couplings between terms where 
l + — I- = constant arc allowed as a consequence of the 
model. Similar to the reduction for a three-mode model 
in the quasi-2D case treated in [l(| , this allows us to par- 
tition the eigenvalue problem into 2N + 1 smaller prob- 
lems, by introducing the notation ipjf l = 4>h,i + .i- with 

I 



—N<k = l_ — l + <N and ip = for indices such that 
| Z x + 2l + + k\ > N. In terms of the many-body ampli- 
tudes with index fc, the new eigenvalue problem, with the 
relations we have found for the matrix elements, c.f. ©, 
reads 



Since e_ = e + and C+ = C_, we have ci lt i + j_ = Ci lt i_j + 
and di,d,2, are also symmetric with respect to inter- 
changing l + and Z_. 

The eigenvalue problem (|13[) then becomes completely 
symmetric in l + and Z_, i.e., interchanging them leaves 
the equation unaltered and we can assume them, backed 
up with numerical simulations, to be identical for the 
many-body ground state, up to a global phase <j) in the 
amplitudes, that is 



Influence of higher l± states N = 5000, G 3 = 5000 



l + = Z_ or formally Vh,l± : ^>h,i + , 



We can immediately deduce from (TTJJ that the ground 
state at k = will be (nearly) degenerate due to the de- 
coupling of even and odd values of l±. The problem can 
thus be split further to separate the practically degen- 
erate eigenstates (up to the energy of one particle) that 
belong to only even or only odd occupation numbers l\. 
These substantial size reductions allow to solve the eigen- 
value problem numerically and yields a matrix size of ap- 
proximately N 2 x N 2 . Due to the coupling structure of 
the Hamiltonian, this matrix is very sparse and the total 
number of nonzero entries grows only quadratically with 
the particle number 0(N 2 ). 

We remark that the even-odd degeneracy allows for the 
free choice of a phase parameter 9 in the superposition 
of the degenerate eigenstates [18|, which will depend on 
the preparation of the state, and which we have taken to 
be zero in what follows. 



B. Numerical results 

The ground states for different configurations N, G3 
have been computed numerically by finding local minima 
in the energy curve along the variational parameter A. 

Numerical calculations for particle numbers up to 
TV = 20000 confirm that k > states correspond to 
increasingly higher energies, with the lowest difference 
(k = to k = 1) in energy per particle being ap- 
proximately 0(G^ AA /NiS) and relative energy differences 
1 - E k =x/E k=Q also scaling with 0{l/N). 
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FIG. 1. (color online) The left panel shows the error commit- 
ted by truncation of the eigenvalue problem at a fixed value 
max l± for different values of the particle number JV and inter- 
action strength G with A m i n fixed at the minimal variational 
energy configuration. Curves with same line style and sym- 
bol correspond to same N or G, respectively. The right plot 
visualizes the locations of amplitudes IV'fi.i+l 2 larger than a 
given threshold. Note the scaling on the h axis which is of 
order O(N), whereas the horizontal axis ends at l± — 7. 



Fixing k at zero, and hence l + — I—, the occupation of 
the circulating orbitals stays below (j ±> o) IV'ii^il 2 < 
3%, with all significant amplitudes located at l± ~ 0(1), 
whereas the occupation of the radially symmetric orbital 
at l± = is scaling with li ~ O(N), cf. Fig.[TJ Including 
as few as eight circular states is sufficient to reach ma- 
chine accuracy for all reasonable configurations TV, G3. 

For large values of G3 a second shallow minimum ap- 
pears in the energy landscape, analogous to the quasi- 
1D and quasi-2D cases [I(J, and the condensate starts 
to fragment. We have determined numerically that the 
onset of fragmentation is determined by a critical value 
of the interaction parameter, (Gs) c , which depends on 
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FIG. 2. (color online) Degree of fragmentation (circles) of the 
ground state for varying N at fixed G3 = 5000 for nonfrag- 
mented (blue empty circles) and fragmented (red filled circles) 
states. Small N effects include a nonvanishing fragmentation 
which quickly approaches 0. After passing a critical N, a new 
minimum appears and asymptotes its maximum for moderate 
values of N. The dashed-dotted black line shows the degree 
of fragmentation at T = 0.16, computed in the large N limit 
(|25l) . Vertical bars (yellow) indicate the sensitivity of the de- 
gree of fragmentation when we allow for an energy variation 
away from the local minimum via A, up to the energy barrier 
separating the two minima. 



FIG. 3. (color online) Dependence of degree of fragmenta- 
tion (circles) on G3 for fixed N = 20000. The single exist- 
ing minimum (blue empty circles) for subcritical G3 shows no 
fragmentation, whereas the one created after a critical G3 has 
been passed (red filled circles) soon asymptotes to its limiting 
value J- — 0.19 indicated by the dashed line. The solid line 
corresponds to the continuum limit results computed via the 
minimization of (|25[) . 



IV. THE CONTINUUM LIMIT OF THE 
EFFECTIVE TWO-MODE MODEL 



particle number. For small N ~ 1000, the critical inter- 
action strength is (G^) c « 5600, a value that decreases 
quickly to its asymptotic value {Gs) c = 2480 ± 10 for 
N = 50000. For large values of G3, the fragmented local 
minimum becomes a global one and the non-fragmented 
minimum becomes very shallow. 

A finite particle number effect on the fragmentation 
can be observed by varying N for given G3 and is il- 
lustrated in Figure [5] For small particle numbers, cor- 
responding to a subcritical Ga(iV), only one minimum 
exists and its degree of fragmentation, defined by J- — 
1 — |Ai — A2I/-/V, where Xi are the (macroscopic) eigenval- 
ues of the single-particle density matrix rapidly ap- 
proaches zero when N is increased. Once we have passed 
the critical value for N, fragmentation sets in, with the 
appearance of a new local energy minimum at a smaller 
extension A, and quickly approaches a limit which de- 
pends on the chosen value of G3. 

The dependence of the fragmentation on the interac- 
tion strength G3 is depicted in Fig. [3] Note the onset 
of fragmentation after passing the critical (G3) C (N) — 
O(10 3 ). 

The considerations above on the low occupancy of cir- 
culating states, cf. Fig. [1] allow to simplify the problem 
by neglecting the small l± > contribution. This yields 
a two-mode model, for which we apply a continuum limit 
[ij} in the following section. 



A. Derivation of the Schrodinger equation for the 
mode population 

From the full numerical analysis, we are led to con- 
clude that l± is approximately unpopulated and af- 
ter deletion of the l± > terms, the eigenvalue problem 
(fT4")) reduces to 

+ ^(M^+a.o + - 2,0)^_ 2j0 . (15) 

The approximations (omitting the subscript 1) d\{l) 
di = -(I - N/2) 2 + N 2 /A and d x (l + 2) w d x {l) yield 
di(l) (ipi+2 — 2ipi + ipi-2) ~ 4di(Z)<9; 2 , which we use to 
write (|15[) as an ordinary differential equation 

4^dx(i)3V(0 + c(l) + 2^d 1 {l)m = E^(l)- (16) 

We employ a change of variables, neglecting from here 
on 0(1/N) terms, t = l- N/2, and obtain 

2E 1 [-t 2 + N 2 /4]d^(t)+ 
[E^t 2 + N 2 /A) + c(t + N/2)] = £7*(t). (17) 

We identify the above equation with the Schrodinger 
equation of the harmonic oscillator via 

-2^* 2 *(*) + Q mw - & f + e ©) *(*) = B *(*)> 

(18) 
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with the parameters 
1 

771 



-Ex{C +C 1 -D 1 - 2E 1 )(N 2 - At 
{e a -e 1 ) + {C -C 1 ){N-l)/2 



2 1 

fix/' 



(19) 
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/V 2 



CJV/2 



-muj 2 e 2 , 



and an energy shift ee independent of t. Note that we 
have fixed the contribution from the mode-exchange (oc 
Ei) at t = tfi x , that is taking <ii(ig x ) instead of e?i(i), in 
front of the derivative. With the additional substitutions 
e = XNC and the scalings t = T ■ N/2, t fix = T fix TV/2 



^AT in terms 



and -Ei = \C$, we then obtain © = — 1 
of the ratio of single-particle energy to interaction energy 
units X = e /(TVG ). 

The ground state of (fTS)) can be solved for analyti- 
cally when the absolute value of the Fock-state ampli- 
tudes \ip(t) \ is considered as a continuous variable [9T.ll7j|. 
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2a 2 



(20) 



We get the effective oscillator length of the "harmonic 
oscillator" (note that T € [—1,1]) as 



a 2 = 



= TV 



3(1 -TD- 



(21) 



The single-particle to interaction energy units ratio X is 
then calculated to be 



X = 



eo 
NCn 



3 _L_ 

4ft 



(A + A 5 



(22) 



Finally, the total energy in the continuum limit for the 
reduced model is given by 
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ee 



Ei 



TV 2 



C/V/2 



- -muj 2 e 2 . 



(23) 



We note that, to this order, the dependence of er on Tfi x 
does not enter the continuum energy; we finally obtain, 
to first order in TV, 



E 
TVCo 



TV 
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and the "shift" becomes & 
u ° 1va?> 
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(24) 
Then, with 

3A 3 ' 
(25) 



which represents the continuum expression for the energy 
as a function of A, with the sole parameter G3. 



B. Large coupling limit 

The minimization problem OE/dA = in the limit of 
can be solved for real values of G3 when A > 4.686, 
and we get for the minimum 
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G 3 = H(A 5 - A) + ^v / 193A 2 - 482A 6 + A 1 
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For the relative energy difference to a Fock state with 
all particles occupying the I = state (the radial ground 
state), we get 



AE 



Fock 



Eco 



Nuj 

with the Fock state energy 



0.00013 G 



2/5 



(27) 



(28) 



At its minimum, Aj> ockmin 



G3. Note also that in 

2/3 



quasi- ID, we had AE/{Nuj) ~ 0.02G7 J and m quasi- 

2D, AB/fiVwi) - 0.002G2 72 @. Hence, with increas- 
ing dimension, both the prefactor as well as the scal- 
ing of the energy difference to a single condensate de- 
crease. For a second, fragmented minimum to exist we 
need G 3 > 2436.13, and then have AE/(Nuj) < 0.0015. 

The continuum limit is valid around the expansion 
point ifi x , which we put equal to the shift, 2fi x = 
&/{N/2) = i - +fX. For the single-particle to 

interaction-energy ratio, we have X = | (f) 1 ^ 5 G^ 4 ^ 5 + 
jt, which asymptotes to X = jt m 0.64 and is close to 
this value already for the critical {G^c = 2436.13. Then, 
the asymptotic shift is evaluated to ©/(TV/2) w —0.810, 

giving cr 2 = N^Z(1 - T| x ) w 0.48TV. We can now assess 

the validity of the continuum approach by measuring the 
occupation it assigns to negative (unphysical) l\. With 
increasing TV, the width of the wavefunction ([20]) gets 
smaller, and the density at negative l\ goes to zero as 
i(l — erf [0.13729-\/TV])- The rapid convergence is illus- 
trated for TV = 1000, when f~* /2 \ijj{t)\ 2 dt w 4.1 • lO" 10 . 



C. Degree of fragmentation 

The degree of fragmentation in the continuum limit 
reads 



N 



TV 



■ sint 



1 - 



26 2 



TV 2 



6 2 



(29) 



Here, we assume that the two degenerate many-body 
states of the two- mode problem 0, have equal weight 
in the ground state, and 9 is their relative phase (l8j . 
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The maximal degree of fragmentation (that is when 
6 = as assumed in our numerical computations above) 
becomes 

Within the validity of the continuum approximation, in 
the limit of large coupling, the fragmentation reaches ap- 
proximately 19%, and is hence significantly lower than 
in either quasi-lD (80%) and quasi-2D (33%) trapping 
geometries. The power law of the asymptotics here is 
4/5, while in quasi-lD trapping it has been 4/3 and in 
quasi-2D unity This implies that the coupling de- 
pendence of the degree of fragmentation becomes weaker 
with increasing dimension. 

Finally, we conclude from the comparison with the nu- 
merical data shown in Figs. [2] and G2 that the agreement 
of two-mode continuum limit and numerics is excellent 
for sufficiently large values of G3 and N. 

V. CONCLUSION 

A dimensionless measure, Gd, of the relative impor- 
tance of total interaction and potential energies, which 
ultimately determines the first-order coherence proper- 
ties of a trapped system, can be constructed from the 
three-dimensional coupling constant g and the relevant 
trapping lengths in quasi-lD, quasi-2D and proper three- 
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